/***************************************************************************
 *      particle.cpp
 *		
 *      Copyright 2012 Artem Korzhimanov <korzhimanov.artem@gmail.com>
 *      
 *      This program is free software; you can redistribute it and/or modify
 *      it under the terms of the GNU General Public License as published by
 *      the Free Software Foundation; either version 2 of the License, or
 *      (at your option) any later version.
 *      
 *      This program is distributed in the hope that it will be useful,
 *      but WITHOUT ANY WARRANTY; without even the implied warranty of
 *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *      GNU General Public License for more details.
 *      
 *      You should have received a copy of the GNU General Public License
 *      along with this program; if not, write to the Free Software
 *      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 *      MA 02110-1301, USA.
 ***************************************************************************/

#include "particle.h"
#include <cmath>

Particle::Particle()
{
	r = new double[3];
	p = new double[3];
	m = 1.;
	q = 1;
}

Particle::Particle(double* position, double* momentum, double mass, double charge)
{
	r = new double[3];
	p = new double[3];
	for (int i = 0; i < 3; i++)
	{
		r[i] = position[i];
		p[i] = momentum[i];
	}
	m = mass;
	q = int(charge);
}

Particle::~Particle()
{
	delete[] r;
	delete[] p;
}

void Particle::SetPosition(double *position)
{
	for (int i = 0; i < 3; i++)
		r[i] = position[i];
}

void Particle::SetMomentum(double *momentum)
{
	for (int i = 0; i < 3; i++)
		p[i] = momentum[i];
}

double Particle::GetGamma()
{
	return sqrt(1. + p[0]*p[0] + p[1]*p[1] + p[2]*p[2]);
}

void Particle::MakeStep(double dt, double *E, double *B)
{
	double F;
	double gamma = GetGamma();
	for (int i = 0; i < 3; i++)
	{
		F = E[i];
		switch (i)
		{
			case 0 : F += (p[1]*B[2] - p[2]*B[1]) / gamma; break;
			case 1 : F += (p[2]*B[0] - p[0]*B[2]) / gamma; break;
			case 2 : F += (p[0]*B[1] - p[1]*B[0]) / gamma; break;
		}
		F *= -q/m;
		r[i] += (p[i] / (gamma*m) + .5 * F * dt) * dt;
		p[i] += F * dt;
	}
}
